#Cleaned Version Of Code - Produces all Attention Analyses
rm(list=ls())
## setwd("")
library(xts)
library(stringr)
library(lubridate)
library(lmtest)
library(nlme)
library(ggplot2)
library(scales)
library(sandwich)

#Data loading
## load("AttentionConsistent.RData")
combined <- read.csv("AttentionConsistent.csv")
combined$Date <- with(combined, as.Date(paste(year(Date), month(Date), "01", sep = "-")))
## combined[,"Date"]<-as.Date(combined[,"Date"])
combined<-xts(combined[,c("Product","IsFopo","Duration","indic","words")],order.by=combined[,"Date"])
finalAttention<-apply.monthly(combined[,c("Product","IsFopo","Duration","indic","words")],colSums,na.rm=T)
#remove the zero durations - these are nonsense
finalAttention<-finalAttention[finalAttention[,"Duration"]>0,]
#A handful of garbage dates fall outside the window
finalAttention<-finalAttention["1933-01-01/1993-01-31"]

#Put in an NA so the plot doesn't join from Truman to Johnson
newentry<-xts(rbind(c(NA,NA,NA,NA,NA)),order.by=as.Date("1962-03-01"))

#Build augmented data frame for easy analysis
augmented<-cbind(finalAttention,month(finalAttention),year(finalAttention))
colnames(augmented)[6:7]<-c("Month","Year")
iselection<-augmented[,"Year"]%%4==0
midterm<-augmented[,"Year"]%%4==2
augmented<-cbind(augmented,iselection,midterm)
colnames(augmented)[8:9]<-c("iselection","ismidterm")
proximate<-augmented[,"Month"]%in%c(6,7,8,9,10,11)&augmented[,"iselection"]
augmented<-cbind(augmented,proximate)
colnames(augmented)[10]<-"proximate"
premidterm<-augmented[,"Month"]%in%c(6,7,8,9,10,11)&augmented[,"ismidterm"]
augmented<-cbind(augmented,premidterm)
colnames(augmented)[11]<-"premidterm"
inWindow<-augmented[,"Month"]%in%c(6,7,8,9,10,11)
augmented<-cbind(augmented,inWindow)
colnames(augmented)[12]<-"inWindow"
augmented<-cbind(augmented,NA)
#Fixed effects are assembled at the President-term level - i.e., switch to a new FE
#When President changes (incl. via death) or after a subsequent innauguration
colnames(augmented)[13]<-"Term"
augmented["1933-03-01/1937-01-31","Term"]<-1
augmented["1937-02-01/1941-01-31","Term"]<-2
augmented["1941-02-01/1945-03-31","Term"]<-3
augmented["1945-04-01/1949-01-31","Term"]<-4
augmented["1949-02-01/1953-01-31","Term"]<-5
augmented["1953-02-01/1957-01-31","Term"]<-6
augmented["1957-02-01/1961-01-31","Term"]<-7
augmented["1961-02-01/1963-11-31","Term"]<-8
augmented["1963-12-01/1965-01-31","Term"]<-9
augmented["1965-02-01/1969-01-31","Term"]<-10
augmented["1969-02-01/1973-01-31","Term"]<-11
augmented["1973-02-01/1974-08-01","Term"]<-12
augmented["1974-08-01/1977-01-31","Term"]<-13
augmented["1977-02-01/1981-01-31","Term"]<-14
augmented["1981-02-01/1985-01-31","Term"]<-15
augmented["1985-02-01/1989-01-31","Term"]<-16
augmented["1989-02-01/1993-01-31","Term"]<-17

## remove noisy diaries
finalAttention2<-finalAttention[finalAttention[,"Duration"]>1000,]

## pdf("FopoProportion_etc.pdf", width = 10, height = 3)
## Figures 3-4
g <- ggplot(finalAttention2) +
  geom_line(aes(x=index(finalAttention2), y=Duration/60))
g <- g + annotate(
  "rect",
  xmin=as.Date(paste(seq(1936, 1988, 4), "-01-01", sep = "")),
  xmax=as.Date(paste(seq(1936, 1988, 4), "-12-31", sep = "")),
  ymin=-Inf, ymax=+Inf, fill=alpha('blue', 1/20)
)
g <- g + xlab("Date") + ylab("Total Monthly Duration (hours)")
g <- g + theme(panel.background = element_rect(fill = 'white', color = 'white'))
plot(g)
##
g <- ggplot(finalAttention) +
  geom_line(aes(x=index(finalAttention), y=indic/words))
g <- g + geom_hline(yintercept=1,col="red",lty=2)
g <- g + annotate(
  "rect",
  xmin=as.Date(paste(seq(1936, 1988, 4), "-01-01", sep = "")),
  xmax=as.Date(paste(seq(1936, 1988, 4), "-12-31", sep = "")),
  ymin=-Inf, ymax=+Inf, fill=alpha('blue', 1/20)
)
g <- g + ylim(c(0, 1))
g <- g + xlab("Date") + ylab("Prop. of Words in Google Corpus")
g <- g + theme(panel.background = element_rect(fill = 'white', color = 'white'))
plot(g)
##
g <- ggplot(finalAttention2) +
  geom_line(aes(x=index(finalAttention2), y=Product/Duration))
g <- g + annotate(
  "rect",
  xmin=as.Date(paste(seq(1936, 1988, 4), "-01-01", sep = "")),
  xmax=as.Date(paste(seq(1936, 1988, 4), "-12-31", sep = "")),
  ymin=-Inf, ymax=+Inf, fill=alpha('blue', 1/20)
)
g <- g + ylim(c(0, 0.5))
g <- g + xlab("Date") + ylab("Foreign Policy Duration / Total Min.")
g <- g + theme(panel.background = element_rect(fill = 'white', color = 'white'))
plot(g)
## dev.off()

finalAttention2$year <- year(index(finalAttention2))
finalAttention2$month <- month(index(finalAttention2))
finalAttention2 <- cbind(finalAttention2, NA)
colnames(finalAttention2)[8]<-"Term"
finalAttention2["1933-03-01/1937-01-31","Term"]<-1
finalAttention2["1937-02-01/1941-01-31","Term"]<-2
finalAttention2["1941-02-01/1945-03-31","Term"]<-3
finalAttention2["1945-04-01/1949-01-31","Term"]<-4
finalAttention2["1949-02-01/1953-01-31","Term"]<-5
finalAttention2["1953-02-01/1957-01-31","Term"]<-6
finalAttention2["1957-02-01/1961-01-31","Term"]<-7
finalAttention2["1961-02-01/1963-11-31","Term"]<-8
finalAttention2["1963-12-01/1965-01-31","Term"]<-9
finalAttention2["1965-02-01/1969-01-31","Term"]<-10
finalAttention2["1969-02-01/1973-01-31","Term"]<-11
finalAttention2["1973-02-01/1974-08-01","Term"]<-12
finalAttention2["1974-08-01/1977-01-31","Term"]<-13
finalAttention2["1977-02-01/1981-01-31","Term"]<-14
finalAttention2["1981-02-01/1985-01-31","Term"]<-15
finalAttention2["1985-02-01/1989-01-31","Term"]<-16
finalAttention2["1989-02-01/1993-01-31","Term"]<-17
##
presidentHours <- aggregate(Product ~ Term, finalAttention2, FUN = mean)
names(presidentHours)[2] <- "AvgProduct"
finalAttention3 <- merge(data.frame(finalAttention2), presidentHours, by.x="Term",by.y="Term")
finalAttention3 <- xts(finalAttention3, order.by=as.Date(paste(finalAttention3$year, finalAttention3$month,"01",sep="-")))

finalAttention3 <- merge(
        finalAttention3,
        zoo(,seq(as.Date("1933-03-01"), as.Date("1990-10-31"), by = "month"))
        )

## pdf("FopoHours_fixedEffects.pdf", width=10,height=3)
g <- ggplot(finalAttention2) +
  geom_line(aes(x=index(finalAttention2), y=Product/Duration))
g <- g + annotate(
  "rect",
  xmin=as.Date(paste(seq(1936, 1988, 4), "-01-01", sep = "")),
  xmax=as.Date(paste(seq(1936, 1988, 4), "-12-31", sep = "")),
  ymin=-Inf, ymax=+Inf, fill=alpha('blue', 1/20)
)
g <- g + ylim(c(0, 0.5))
g <- g + xlab("Date") + ylab("Foreign Policy / Total Recorded")
g <- g + theme(panel.background = element_rect(fill = 'white', color = 'white'))
plot(g)
## Figure 5
g <- ggplot(finalAttention3) +
  geom_line(aes(x=index(finalAttention3), y=(Product-AvgProduct)/60))
g <- g + annotate(
  "rect",
  xmin=as.Date(paste(seq(1936, 1988, 4), "-01-01", sep = "")),
  xmax=as.Date(paste(seq(1936, 1988, 4), "-12-31", sep = "")),
  ymin=-Inf, ymax=+Inf, fill=alpha('blue', 1/20)
)
g <- g + xlab("Date") + ylab("Foreign Policy - Term Average")
g <- g + theme(panel.background = element_rect(fill = 'white', color = 'white'))
g <- g + geom_hline(yintercept=1,col="red",lty=2)
plot(g)
g <- ggplot(finalAttention3) +
  geom_line(aes(x=index(finalAttention3), y=(Product/AvgProduct)))
g <- g + annotate(
  "rect",
  xmin=as.Date(paste(seq(1936, 1988, 4), "-01-01", sep = "")),
  xmax=as.Date(paste(seq(1936, 1988, 4), "-12-31", sep = "")),
  ymin=-Inf, ymax=+Inf, fill=alpha('blue', 1/20)
)
g <- g + xlab("Date") + ylab("Foreign Policy / Term Average")
g <- g + theme(panel.background = element_rect(fill = 'white', color = 'white'))
g <- g + geom_hline(yintercept=1,col="red",lty=2)
plot(g)
## dev.off()

## Table 8
reelect<-augmented[,"proximate"]&(!augmented[,"Year"]%in%c(1952,1960,1968,1988))
norelect<-augmented[,"proximate"]&(augmented[,"Year"]%in%c(1952,1960,1968,1988))

model5<-lm(Product/60~reelect+norelect+premidterm+as.character(Month)+as.character(Term),data=augmented)
coeftest(model5,vcov=vcovHAC(model5))
model6<-lm(Product/60~reelect+norelect+premidterm+as.character(Month)+as.character(ceiling(Year/2)),data=augmented)
coeftest(model6,vcov=vcovHAC(model6))


## Table 9
presidentHours <- aggregate(Product ~ Term, augmented, FUN = mean)
names(presidentHours)[2] <- "AvgProduct"
augmented <- merge(data.frame(augmented), presidentHours, by.x="Term",by.y="Term")
augmented <- xts(
    augmented,
    order.by=as.Date(paste(augmented$Year, augmented$Month, "01",sep="-"))
    )

model0<-lm(Product/AvgProduct~scale(Year)*proximate+as.character(Month)+as.character(Term),data=augmented)
coeftest(model0,vcov=vcovHAC(model0))

## Table 10
model00<-lm(Product/60~proximate+as.character(Month)+as.character(Term),data=subset(augmented, index(augmented)>="1946-01-01"))
coeftest(model00,vcov=vcovHAC(model00))
model01<-lm(Product/60~proximate+as.character(Month)+as.character(ceiling(Year/2)),data=subset(augmented, index(augmented)>="1946-01-01"))
coeftest(model01,vcov=vcovHAC(model01))

#BLS Recessions
augmented<-cbind(augmented,F)
colnames(augmented)[ncol(augmented)]<-"Recession"
augmented["1937-05-01/1938-06-01","Recession"]<-TRUE
augmented["1945-02-01/1938-10-01","Recession"]<-TRUE
augmented["1948-11-01/1949-10-01","Recession"]<-TRUE
augmented["1953-07-01/1954-05-01","Recession"]<-TRUE
augmented["1957-08-01/1958-05-01","Recession"]<-TRUE
augmented["1960-04-01/1961-02-01","Recession"]<-TRUE
augmented["1969-12-01/1970-11-01","Recession"]<-TRUE
augmented["1973-12-01/1975-03-01","Recession"]<-TRUE
augmented["1980-01-01/1980-07-01","Recession"]<-TRUE
augmented["1981-07-01/1982-11-01","Recession"]<-TRUE
augmented["1990-07-01/1991-03-01","Recession"]<-TRUE

FEMA <- read.csv("FEMA_AllYears.csv")
FEMA$Date <- as.Date(FEMA$Date, "%m/%d/%Y")

FEMA <- xts(
    FEMA,
    order.by=as.Date(paste(year(FEMA$Date), month(FEMA$Date), "01", sep = "-"))
    )
FEMA <- apply.monthly(FEMA, nrow)
names(FEMA) <- "disasters"


augmented <- merge(
    augmented,
    FEMA
    )

augmented$disasters[is.na(augmented$disasters)&index(augmented)>="1953-01-01"] <- 0

## Table 16
model1<-lm(Product/60~proximate+Recession+as.character(Month)+as.character(Term),data=subset(augmented, index(augmented)<="1990-10-01"))
coeftest(model1,vcov=vcovHAC(model1))

model2<-lm(Product/60~proximate+scale(disasters)+as.character(Month)+as.character(Term),data=subset(augmented, index(augmented)>="1953-01-01" & index(augmented)<="1990-10-01"))
coeftest(model2,vcov=vcovHAC(model2))
